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I .  Introduction 

This  report  describes  three  procedures  for  estimating  the 
position  of  an  object.  Each  procedure  is  based  on  a  model  that 
describes  the  position  of  an  object  in  terms  of  the  position  of  a 
point  on  a  plane  surface. 

The  first  procedure  is  developed  in  Section  II.  The 
procedure  is  a  method  for  combining  into  a  composite  estimate  two 
or  more  point  estimates  of  the  position  of  an  object  that  each 
have  confidence  regions  that  are  bounded  by  an  ellipse.  In  the 
model  on  which  it  is  based,  the  coordinates  of  each  estimate  are 
the  values  of  random  variables  that  have  a  bivariate  normal 
distribution  with  a  mean  vector  whose  elements  are  the  unknown  ' 
coordinates  of  the  object.  A  user  of  the  procedure  must  specify 
the  coordinates  of  each  estimate  and  the  major  and  minor  axes  of 
the  ellipse  that  bounds  its  confidence  region. 

The  second  procedure  is  developed  in  Section  III.  The 
procedure  is  a  method  for  generating  a  position  estimate  from 
bearings  that  are  on  or  from  two  or  more  reference  points  of 
known  position.  In  the  model  on  which  it  is  based,  bearing 
errors  are  independent  normally  distributed  random  variables. 
Based  on  the  model,  the  coordinates  of  a  position  estimate  that 
is  determined  by  the  procedure  are  values  of  random  variables 
that  have  a  bivariate  normal  distribution  with  a  mean  vector 
whose  elements  are  the  unknown  coordinates  of  the  object.  A  user 
of  the  procedure  must  specify  the  standard  deviations  of  the 
bearing  errors.  An  extension  of  the  procedure  is  also  developed 
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in  this  section.  It  is  a  method  for  generating  a  position 
estimate  from  bearings  that  are  on  or  from  two  or  more  reference 
points  for  which  the  position  of  one  or  more  of  the  reference 
points  is  uncertain. 

The  third  procedure  is  described  in  Section  IV.  The 
procedure  is  a  method  for  generating  a  position  estimate  that  is 
determined  from  lines  of  position.  In  the  model  on  which  the 
development  is  based,  lines  of  position  are  straight  lines,  lines 
of  position  chat  arc  determined  by  observation  are  parallel  to 
true  lines  of  position  and  the  algebraic  distance  between  a  line 
of  position  based  on  an  observation  and  a  true  line  of  position 
is  the  value  of  a  normally  distributed  random  variable.  A  user 
of  the  procedure  must  specify  the  standard  deviations  of  these 
random  variables. 

In  the  models  on  which  the  three  procedures  are  based, 
systematic  errors  represent  mean  values  that  are  set  to  zero. 

This  suggests  that  the  procedures  should  be  used  only  were  the 
bias  in  the  measurements  that  determine  bearings  or  lines  of 
position  is  known  or  is  known  to  be  negligible. 
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II.  A  Composite  Position  Estimate 

The  procedure  that  is  developed  in  this  section  is  for 
combining  position  estimates  for  an  object  that  are  generated 
from  data  from  independent  sources.  It  is  based  on  the  following 
model:  In  a  rectangular  coordinate  system,  the  coordinates  of 

each  position  estimate  are  values  of  random  variables  that  are 
determined  by  an  independent  bivariate  normal  distribution  with  a 
covariance  matrix  whose  elements  are  known  and  with  a  mean  vector 
whose  elements  are  the  unknown  coordinates  of  the  object.  The 
above  model  implies  that  the  natural  logarithm  of  the  likelihood 
function  L  for  a  set  of  n  estimates  can  be  expressed  as 
follows:  log  L  =  K  -  (x,  -  x)  Z-1  (xi  -  x)  where  K  is  a 

constant,  x(  is  an  estimate  of  the  mean  vector  with  components 
xi  and  y(  ,  x  is  the  mean  vector  with  components  x  and  y,  the 
unknown  coordinates  of  the  object,  and  Z,  is  the  covariance 
matrix  with  elements  a ,  a\-  and  The  composite  maximum 

likelihood  estimates  x  and  y  of  the  unknown  coordinates  x 
and  y  make  loq  L  a  maximum  which  implies  the  following 
relations: 


After  evaluating  the  derivatives,  they  can  be  written  as: 
A  x  +  B  y  =  D 
B  x  +  C  y  =  E 
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where: 

A  =  Si  a,,  B  =  2j  b(  ,  C  =  c,  ,  D  =  (aixi  +  b^)  , 

E  =  2"  (bixi  +  c,^)  ,  ai  -  o2-  /  di  ,  bi  =  -  aiX-  /  di  , 

c.  =  oz-  /  d.  .  d  =  <j2-  a2-  -  o 

The  solutions  to  the  equations  are: 

X  =  (CD  -  BE)  /  (AC  -  B2) 
y  =  (AE  -  BD)  /  (AC  -  B2) 

Since  the  estimates  x  and  y  are  linear  combinations  of  the 
estimates  xi  and  yf  ,  they  are  determined  by  a  bivariate 
normal  distribution  that  is  defined  by  its  mean  vector  with 
components  x  and  y  and  its  covariance  matrix  with  elements 

o\,  a-  and  a-.  With  E  the  expected  value  operator: 

x  =  E {  (CD  -  BE) /(AC  -  B2)  } 

y  =  E  {  (AE  -  BD)  /  (AC  -  B2)  } 

<7*  =  E{(CD  -  BE)  -  E  (CD  -  BE))2/ (AC  -  B2)2 

=  (C2(F+I+2L)  -  2CB (H+K+M)  +  B2  (G+J+2N)  }/ (AC  -  B2)2 

a]  =  E  {  (AE  -  BD)  -  E  (AE  -  BD))2/(AC  -  B2)2 

=  (A2  (G+J+2N)  -  2AB  (H+K+M)  +  B2  (F+I+2L)  )/ (AC  -  B2)2 

ai9  =  E{[(cd-be)-E(cd-be)][(ae-bd)-E(ae-bd)])/(ac  -  b2)2 

=  (  (AC+B2)  (H+K+M)  -CB(F+I+2L)  -BA(G+J+2N)  )/ (AC  -  B2)2 
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where : 


F  -  2"  a 2a2-  ,  G  —  2"  b24  ,  H  -  ,  I  -  b^c-  , 


n  u.2^2 


j  -  j;  -  K  -  s;  b,c,4,  ,  L  =  s;  a(b,a)W  ,  M  -  2;  (a,c,+b5). 


ixy 


N  =  2;  b jCial5- 


Without  loss  of  generality,  the  origin  of  the  coordinate 
system  can  be  the  point  determined  by  the  composite  estimate,  its 
positive  y-axis  can  be  north  and  its  positive  x-axis  east.  In 
Section  II  of  this  report,  it  is  shown  that  the  axes  of  the 
ellipse  that  bounds  a  minimum  area  confidence  region  for  a 
composite  estimate  are  coincident  with  a  primed  coordinate  system 
which  is  transformed  to  this  coordinate  system  by  a  coordinate 
axes  rotation  equal  to  7  where  7  is  defined  by  the  relation: 
tan  27  =  2 a*-  /(a?  -  a?)  . 

It  can  be  shown  that  a  minimum  area  confidence  region  for  a 
two  dimensional  estimator  whose  components  have  a  bivariate 
normal  distribution  is  bounded  by  an  ellipse  whose  nujor  axes  are 

t  1  /? 

determined  by  2k o*.  and  2ka-.  where  k  =  [-2  log(l  -  p)  ] 
and  x'  and  y*  refer  to  a  coordinate  system  whose  axes  are 
parallel  to  the  axes  of  the  bounding  ellipse.  In  this  system: 


4  = 

4 

cos2 

7  - 

2 

xy 

sin 

7 

cos 

7 

+ 

4  = 

4 

sin2 

7  + 

2 

a— 

xy 

sin 

1 

cos 

1 

+ 

1  .  2 

o ^  sm^  7 

2  2 

a -  cos  7 
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The  model  that  is  described  above  can  be  used  to  find  a 
composite  position  estimate  from  two  or  more  position  estimates 
for  an  object.  To  show  this,  first  assume  that  the  estimates 
satisfy  the  conditions  that  are  required  by  the  model.  Then,  let 
5,  be  the  direction  of  the  major  axis  of  the  ellipse  that  bounds 
the  confidence  region  that  is  associated  with  the  i  th  position 
estimate.  Next,  orient  the  primed  coordinate  system  defined 
above  so  that  its  x'-axis  is  coincident  with  the  major  axis  of 
the  ellipse.  Then,  the  elements  of  the  covariance  matrix  of  the 
bivariate  normal  distribution  that  determines  the  estimate  are 
given  by  the  following  relations: 


2  2  -2 

sln 

+ 

2 

COS 

(J-  =  COS2 

ly  lx 

6i  + 

2 

sm 

a =  (a2-.  - 

lxy  '  ix 

sin 

& ,  cos  6 

As  an  example  of  the  use  of  the  procedure,  suppose  that  each 
bounding  ellipse  were  a  circle.  In  this  case,  the  coordinates  of 
each  position  estimate  would  be  determined  by  a  circular  normal 
distribution  and,  for  i  =  1  to  n  ,  the  elements  of  the  covariance 
matrices  of  the  circular  normal  distributions  would  satisfy  the 
relations:  a-  =  a i-  =  ai  and  o-  =  0.  The  composite  estimate 

would  be: 


[2? 

(*i 

/  a,)] 

/ 

[2"  (l  / 

®i)] 

[2? 

<y. 

/ 

/ 

t  (2J  1  / 

a,)] 
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And  the  elements  of  the  covariance  matrix  of  the  distribution 
that  determine  the  coordinates  of  the  estimate  would  be: 

o\  =  n  /  [ZJ  (1  /  ai )  ] 2  ,  a?  =  n  /[(ZJ  (1  /  a,)]2  and  ok-  =  0  . 

Consequently,  the  distribution  would  also  be  circular  normal  and 
the  bounding  curves  for  minimum  area  confidence  regions  would  be 
circles . 

A  program  called  COMP  that  is  based  on  the  procedure  that  is 
developed  in  this  section  is  listed  in  Appendix  1.  COMP  can  be 
used  to  generate  composite  position  estimates  and  corresponding 
confidence  regions.  As  an  example  of  its  use,  suppose  three 
independent  position  estimates  are  given  in  Table  1. 


i 

*1 

y  i 

TABLE  1 

MJi 

MIf 

k 

1 

-3.7 

18.1 

59° 

36 

20 

2 

2 

11.8 

00 

4* 

105° 

38 

10 

2 

3 

0 

0 

146° 

50 

24 

2 

In  Table  1,  MJ  is  the  length  of  the  major  axis  and  MI  is  the 
length  of  the  mir. ->r  axis  of  the  ellipse  bounding  a  minimum  area 
confidence  region.  If  the  distance  unit  is  the  nautical  mile, 
the  coordinates  of  the  composite  estimates  of  the  coordinates  of 
the  object  are:  x  =  -2.69  nm  and  y  =  12.41  nm.  For  k  =  2, 
the  lengths  of  the  m-jor  and  minor  axes  of  the  bounding  ellipse 
of  the  minimum  area  confidence  region  are:  MJ  =  17.33  nm  and 
MI  =  8.85  nm  and  the  direction  of  its  major  axis  is  103.77°. 
The  results  are  shown  in  Figure  1. 
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III.  A  Position  Estimate  Based  on  Bearings 

The  procedure  that  is  developed  in  this  section  is  for 
determining  a  position  estimate  for  an  object  from  a  set  of 
observed  bearings  from  or  on  two  or  more  reference  points  of 
known  positions.  In  the  development,  reference  points  are 
referred  to  as  stations.  The  procedure  is  based  on  the  following 
model:  A  set  of  bearings  taken  on  an  object  from  a  station  are 
the  values  of  independent  normally  distributed  random  variables 
with  known  standard  deviations  but  with  means  that  are  equal  to 
the  unknown  true  bearings.  The  model  implies  that  the  natural 
logarithm  of  the  likelihood  function  L  for  a  set  of  n 
bearings  taken  on  an  object  can  be  expressed  as  follows: 
log  L  =  K  -  1/2  (0,  -  <f>.)2  /  e^  where  K  is  a  constant,  8- 

is  an  observed  bearing,  <f>s  is  the  unknown  true  bearing  and  e, 
is  the  known  standard  deviation.  For  i  =  1  to  n  ,  the  maximum 
likelihood  estimates  0,  are  the  solutions  of  the  n  equations: 

agog  L)  „  =  0 

30.  0,=0 , 

subject  to  the  constraint  that  the  bearing  lines  determined  by 
the  estimates  must  all  pass  through  a  common  point.  Because  of 
this  constraint,  the  number  of  independent  estimates  is  reduced 
from  n  to  2.  This  is  equivalent  to  stating  that  any  two  of 
the  estimates  determine  a  common  point.  The  common  point  is  the 
maximum  likelihood  estimate  for  the  position  of  the  object. 
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In  order  to  impose  the  constraint  on  the  estimates,  the 
likelihood  function  will  be  reformulated  in  a  rectangular 
coordinate  system.  To  do  this,  first,  consider  the  following 
definitions: 

u,  =  ri  -  <t>s)  ,  vi  =  r,  (0,  -  0S)  and  wi  =  r,  (0,  -  /J, ) 

where  i  =  1  to  n  is  the  observation  number  and  the  observations 
are  from  two  or  more  stations.  In  these  definitions,  is  the 

bearing  of  an  initial  estimate  of  the  position  of  an  object  and 
ri  is  its  range  from  the  station  associated  with  the  i  th 
observation.  It  is  also  the  radius  of  a  circle  that  is  centered 
on  the  station  and  passes  through  the  initial  estimate  as  shown 
in  Figure  2.  With  0 ,  ,  0,  and  0 t  in  radians,  ui  ,  vi  and  wi 
are  the  length  of  arcs  on  this  circle.  And,  ut  =  w(  -  v(  .  In 
this  relation,  vi  involves  the  unknown  true  bearing  0i  ,  but  wi 
is  known. 

In  order  to  express  vi  in  a  convenient  form,  consider  a 
rectangular  coordinate  system  whose  origin  is  at  the  initial 
estimate,  whose  positive  y-axis  is  north  and  whose  positive 
x-axis  is  east.  With  x  and  y  the  coordinates  of  the  object, 
to  first  order:  v1  =  x  cos  jSi  -  y  sin  .  This  implies: 

=  0*  +  (x  cos  -  y  sin  0,)  /  r, 

And,  to  first  order,  this  expresses  the  constraint  on  the  0i 
and  on  their  maximum  likelihood  estimates  and  it  implies  that  the 
bearing  lines  determined  by  0]  ,  0t  and  $ i  are  parallel. 
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Figure  2.  In  the  program  that  implements  the  procedure,  the 
origin  of  the  xy-coordinate  system  is  coincident  with  the  common 
reference  point  rather  than  with  the  initial  estimate  and  the 
coordinates  of  the  initial  estimate  are  x*  and  y*.  The 
program  is  discussed  later  in  the  section. 
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In  order  to  determine  the  estimates  <f>i  to  first  order,  use 
the  above  relation  to  replace  by  the  known  values  Ps  and 

the  two  unknown  values  x  and  y.  Next,  determine  the  maximum 
likelihood  estimates  x  and  y  for  the  unknown  coordinates  x 
and  y  which,  to  first  order,  are  the  unknown  coordinates  of  the 
object.  That  is,  find  the  solutions  to  the  equations  determined 
by: 

=  0  and  d (loa  L)  =  0 

x=x  dy  x=x 

y=y  y=y 

Evaluating  the  derivatives,  the  equations  are: 


yn 

[  (wi 

X  cos  P- 

+  y 

sin 

*i> 

cos  p.] 

o 

ll 

to" 

\ 

[  (w, 

X  COS 

+  y 

sin 

sin 

/  o  1  =  0 

where  a,  =  xi  e(.  These  equations  can  also  be  written  as: 

A  x  +  B  y  =  D 
B  x  +  C  y  =  E 

and  their  solutions  as: 

X  =  (CD  -  BE)  /  (AC  -  B2) 

y  =  (AE  -  BD)  /  (AC  -  B2) 

where 

A  =  Zj  ( cos2  0()  /  a2  ,  B  =  Zj  (cos  ps  sin  ps)  /  a 2  , 

C  =  Zj  (sin2  ps)  /  a\  ,  D  =  Zj  (wi  cos  )  /  o\ 

E  =  Z"  ( Wj  sin  pi )  /  a]  . 
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If  the  distribution  of  the  random  variables  that  determine  x 
and  y  have  a  known  distribution,  then,  in  principle,  a 
confidence  region  for  an  estimate  can  be  constructed.  With  this 
in  mind,  the  distribution  of  the  estimators  that  determine  x 
and  y  will  be  determined  next.  If  X  and  Y  represent  the 

estimators,  that  is,  the  random  variables  that  determine  x  and 
y  ,  then , 

X  =  (Wi  /  a])  (B  sin  -  C  cos  ps )  /  (B2  -  AC) 

Y  =  (Wf  /  a])  (A  sin  j  ~  B  cos  0f)  /  (B2  -  AC) 

where  Wi  =  ri  (01  -  /?,)  .  Since  the  01  are  normally  distributed 
random  variables,  the  Wi  are  also.  And,  since  the  estimators  X 
and  Y  are  linear  combinations  of  the  W,  ,  they  have  a 
bivariate  normal  distribution. 

If  0,  =  ,  for  i  =  1  to  n,  that  is,  if  the  initial 

estimate  is  the  position  of  an  object  then,  with  E  again  the 

expected  value  operator,  E(W^)  =  0  ,  since  Wi  =  r,  (0i  -  y?i )  .  In 
this  case: 

a\  =  (l/a2i)  (B  sin  -  C  cos  Px )2/(B2  -  AC)2 

0?  =  Zj  (l/0Zi)(A  sin  0 -  B  cos  p^)2/(B2  -  AC)2 

0j-  =  (1/02,)  (B  sin  Ps  -  C  cos  (A  sin  0,  -  B  cos  ^)/(B2  -  AC)2 

Using  the  definitions  for  A  ,  B  ,  and  C  ,  that  are  given 
above,  these  equations  can  be  written  as: 
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a\  =  C  /  (AC  -  B2) 
o\  =  A  /  (AC  -  B2) 

Ory  =  B  /  (AC  -  B2) 

To  the  degree  that  the  initial  estimate  is  near  the  position  of 
an  object,  these  equations  define  approximations  for  the  elements 
of  the  covariance  matrix  of  the  estimators  X  and  Y. 

Since,  in  general,  will  not  be  zero,  the  axes  of  the 

ellipse  that  bounds  a  minimum  area  confidence  region  for  x  and 
y  will  not  be  parallel  to  the  xy-coordinate  axes.  To  determine 
the  orientation,  translate  the  xy-coordinate  system  so  that  its 
origin  is  at  the  point  determined  by  x  and  y  .  Then,  rotate 
this  coordinate  system  until  its  axes  coincide  with  the  axes  of 
the  bounding  ellipse.  Next,  represent  coordinates  in  this  system 
by  x'  and  y'.  If  the  rotation  angle  is  y  ,  then  the 
coordinates  of  a  point  in  the  two  systems  are  related  by  the 
following  general  transformation  equations: 

x'  =  x  cos  y  -  y  sin  y 
y'  =  x  sin  y  -  y  cos  y 

These  transformation  equations  can  be  used  to  determine  the 
elements  of  the  covariance  matrix  of  the  distribution  in  the 
primed  coordinate  system  from  the  covariance  elements  of  the 
distribution  matrix  in  an  unprimed  coordinate  system,  by  using 
the  following  equations: 
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Since  X'  and  Y'  are  independent  random  variables  when  the  axes 
of  the  coordinate  system  are  parallel  to  the  axes  of  the  bounding 
ellipse,  o^-,  =  0  which  implies  that  7  satisfies  the  equation: 

tan  27  =  2 a--  /(a?  -  a\)  . 

If  the  position  of  an  object  were  known,  and  the  estimated 
position  were  at  the  true  position,  then  the  probability  p  that 
an  estimated  position  would  be  in  a  region  bounded  by  an  ellipse 
with  axes  2kas.  and  2k a-,  that  is  centered  on  the  true  position 
is  given  by  p  =  1  -  exp(-k2/2)  .  This  result  can  be  obtained  by 
integrating  over  the  region  bounded  by  the  ellipse  the  bivariate 
normal  density  that  determines  X'  and  Y'.  Note  that  k  is  a 
measure  of  the  size  of  the  region  bounded  by  the  ellipse  and,  for 
a  probability  of  containment  p  ,  the  corresponding  value  of  k 
is  determined  by  k  =  [-2  log(l  -  p) )  ] 1/2  .  These  relations 
provide  a  basis  for  determining  confidence  regions  for  the 
estimates  x  and  y  of  the  position  of  an  object  that  are  given 
above . 

A  program  called  PEST  that  is  listed  in  Appendix  2  can  be 
used  to  generate  position  estimates  and  corresponding  confidence 
regions  from  bearing  observations  on  or  from  two  or  more 
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stations.  For  a  bearing  on  a  station,  8  is  the  reciprocal  of 
the  observed  bearing.  In  the  program,  the  initial  estimate  is 
the  point  of  intersection  of  the  two  bearing  lines  that  are 
determined  by  the  first  two  entered  bearings  with  each  from  a 
different  station.  After  the  first  two  bearings  have  been 
entered,  the  program  determines  the  coordinates  x*  and  y*  of 
the  initial  estimate  by  solving  the  following  equations: 

x*  sin  ( 82  -  8X)  =  [Pj  sin  (q1  -  8^)  ]  sin  8Z 

-  [pz  sin  (a2  -  8Z)]  sin  8 j 

y*  sin  (8Z  -  8 j)  =  [Pj  sin  (Oj  -  8 ]  cos  8Z 

-  [p2  sin  (a2  -  8Z)  ]  cos  8X 

where  o1  and  pl  are  the  bearing  and  range  of  the  first 
station  and  a2  and  p2  are  the  bearing  and  range  of  the  second 
station  from  a  common  reference  point.  This  method  of  finding 
the  initial  estimate  is  based  on  one  that  is  used  in  Reference  3. 
The  point  estimates  a  and  p  of  the  bearing  and  range  of  an 
object  from  a  common  reference  point  are  shown  in  Figure  3. 

As  an  example,  consider  the  data  in  Table  2. 


TABLE  2 

Station 

Observed 

Station 

Station 

Bearing 

Number 

Bearing 

Bearing 

Range 

Error 

1 

38° 

334° 

13500 

4° 

2 

324° 

50° 

11350 

3° 

3 

3° 

0° 

0 

4  ° 
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North 


East 


Figure  3.  A  confidence  region  bounded  by  an  ellipse  and  a  primed 
coordinate  system  in  which  a-.^.  is  zero.  The  center  of  the 
ellipse  and  the  origin  of  the  primed  and  the  unprimed  coordinate 
systems  are  at  the  estimated  position  of  an  object.  The 
direction  of  the  major  axis  of  the  ellipse  is  6  where 
0°  <  6  <  180'.  The  rotation  angle  from  the  unprimed  to  the 
primed  coordinate  system  is  y  .  The  estimates  a  and  p  are 
the  bearing  and  range  of  the  object  from  a  common  reference 
point . 
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The  order  of  data  entry  in  the  program  is  determined  by  the 
station  number.  Note  that  the  reference  point  is  at  the  location 
of  Station  3.  If  the  distance  ir  in  meters,  then,  for  the  three 
observations,  the  estimates  generated  by  the  program  for  the 
bearing  and  range  of  the  object  from  Station  3  are:  a  =  0.15° 
and  p  =  19554  meters.  For  k  =  2,  the  lengths  of  the  major  and 
minor  axes  of  the  ellipse  bounding  the  minimum  area  confidence 
region  are:  MJ  =  3426  meters  and  MI  =  2260  meters  and  the 
direction  of  its  major  axis  is  12.48°. 

An  extension  of  the  above  procedure  is  developed  next.  The 
extension  accounts  for  cases  where  not  all  of  the  positions  of 
the  reference  stations  can  be  considered  to  be  known.  In  these 
cases,  if  the  station  position  were  to  become  known,  then,  for  a 
bearing  line  that  is  determined  by  an  observation,  the  bearing 
line  that  is  drawn  from  an  assumed  station  position  will  be 
parallel  to  the  bearing  line  that  is  drawn  from  the  station 
position.  The  procedure  is  based  on  this  relation  and  on  the 
model  with  the  addition  of  the  following  condition:  The 
coordinates  of  an  uncertain  station  position  are  determined  by  a 
distribution  wh':ch  is  bivariate  normal  with  a  covariance  matrix 
whose  elements  are  known  and  a  mean  vector  whose  elements  are  the 
unknown  coordinates  of  the  station.  This  implies  that  for  a 
bearing  determined  by  an  observation,  the  algebraic  distance 
between  the  bearing  line  from  an  assumed  station  position  and  the 
bearing  line  from  the  station  position  is  a  random  variable  S 
that  has  a  normal  distribution  with  mean  zero  and  a  standard 
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deviation  as  that  is  determined  by  the  elements  of  the  known 
covariance  matrix. 

The  geometry  of  the  model  is  shown  in  Figure  4.  There,  8 
is  an  observed  bearing  and  the  common  origin  of  the  coordinate 
systems  is  the  assumed  station  position.  For  the  unprimed 
system,  the  positive  y-axis  direction  is  north  and  the  positive 
x-axis  direction  is  east.  For  the  primed  system,  the  x'-axis 
is  coincident  with  the  major  axes  of  an  ellipse  that  bounds  a 
confidence  region  for  the  assumed  position.  The  direction  of  the 
major  axis  is  S  where  0 °  <  6  <  180°.  For  the  double  primed 
system,  the  y”-axis  is  coincident  with  the  assumed  bearing  line. 

The  effect  of  station  position  uncertainty  is  accounted  for 
by  the  algebraic  distance  s  between  a  bearing  line  from  an 
assumed  station  position  and  the  parallel  bearing  line  from  the 
station  position.  In  effect,  to  first  order,  the  following 
bearing  lines  are  also  parallel:  the  bearing  line  from  a  station 
to  the  initial  estimate,  the  bearing  line  from  the  station  to  an 
object  and  a  bearing  line  determined  by  an  observation.  Hence, 
the  random  variable  Uj  =  ri  (0S  -  0,)  is  the  algebraic  distance 
between  the  observed  bearing  line  and  the  bearing  line  from  the 
station  to  an  object.  And  Ue1  =  ri  (0S  -  ±  S  is  a  random 

variable  that  determines  the  algebraic  distance  between  the 
observed  bearing  line  from  an  uncertain  station  position  and  the 
bearing  line  from  the  station  position  to  the  object.  Since  U, 
and  S  are  independent  random  variables,  the  standard  deviation 
of  uei  is:  aei  =  [a]  +  a*]1/z  . 
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Figure  4.  The  geometry  associated  with  the  determination  of  the 
standard  deviation  as  .  Here,  6  is  an  observed  bearing  and 
the  common  origin  of  the  coordinate  systems  is  the  assumed 
station  position.  The  positive  y-axis  direction  is  north  and 
the  positive  x-axis  direction  is  east.  The  x'-axis  is 
coincident  with  the  major  axes  of  an  ellipse  that  bounds  a 
confidence  region  for  the  assumed  station  position,  <5  is  the 
direction  of  its  major  axis  and  the  y"-axis  is  coincident  with 
the  assumed  location  of  the  observed  bearing  line. 


In  the  extended  procedure,  for  a  station  of  uncertain 
position,  the  standard  deviation  a i  =  riei  is  replaced  by: 

°ei  =  [ff?  + 

where 

a2  =  al .  sin2  fS  -  0)  +  a2.  cos2  (6  -  6) 

based  on  Figure  4. 

In  Figure  4,  note  that  ka2.  and  ka2.  are  the  lengths  of 
the  major  and  minor  axes  of  the  bounding  ellipse  for  some  value 
of  k.  For  an  uncertain  station  position,  ^  the  bearing  and 
rt  the  range  of  the  initial  estimate  from  the  station  position 
are  determined  by  the  assumed  station  position.  These 
relationships  are  illustrated  in  Figure  5. 

In  the  development,  although  the  position  of  one  or  more  of 
the  stations  may  be  unknown,  the  position  of  a  common  reference 
point  is  known.  This  implies  that  the  assumed  station  positions 
and  the  position  of  the  initial  estimate  are  also  known. 


Figure  5.  In  the  development  for  the  extended  procedure,  /?  and 
r  are  the  range  and  bearing  of  the  initial  estimate  from  the 
assumed  station  position. 


IV.  A  Position  Estimate  Based  on  Lines  of  Position 

The  procedure  that  is  described  in  this  section  is  for 
determining  a  position  estimate  for  an  object  from  two  or  more 
lines  of  position  that  can  be  considered  to  be  straight  lines  in 
the  neighborhood  of  the  object.  The  procedure  is  based  on  the 
following  model:  A  line  of  position  is  a  straight  line.  A  line 
of  position  that  is  based  on  an  observation  is  parallel  to  the 
true  line  of  position.  The  algebraic  distance  u  between  a  line 
of  position  that  is  based  on  an  observation  and  the  true  line  of 
position  is  the  value  of  a  normally  distributed  random  variable 
with  zero  mean  and  known  standard  deviation. 

Labeling  lines  of  position  by  the  index  i  ,  the  model 
implies  that  the  algebraic  distance  wi  between  the  i  th  line  of 
position  based  on  an  observation  and  a  parallel  line  through  an 
initial  estimate  of  an  unknown  position  is  related  to  the 
algebraic  distance  ui  between  it  and  the  true  line  of  position 
by:  ui  =  Wj—  vi  where  v1  is  the  algebraic  distance  between  the 

i  th  true  line  of  position  and  the  parallel  line  through  the 
initial  estimate.  Now,  consider  a  rectangular  coordinate  system 
whose  origin  is  at  the  initial  estimate,  whose  positive  y-axis 
is  north  and  whose  positive  x-axis  is  east.  With  x  and  y 
the  coordinates  of  the  unknown  position,  to  first  order: 

Vj  =  x  cos  -  y  sin  ^ 

where  is  the  direction  of  the  i  th  true  line  of  position  and 

O'  <  /Jj  <  180".  By  referring  to  the  development  in  Section  III, 
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note  that  the  estimates  x  and  y  that  are  developed  there  also 
apply  here.  However,  in  this  development  the  standard  deviation 
a  of  the  random  variable  that  determines  the  distance  u  must  be 
known  for  each  of  the  observed  bearing  lines. 

Figure  6  illustrates  the  geometry  associated  with  the 
development. 
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Figure  6.  The  geometry  for  three  lines  of  position.  The 
direction  of  the  lines  are  ,  jS2  and  jS3.  With  the  initial 
estimate  of  the  unknown  position  at  the  intersection  of  line  1 
and  line  2  ,  both  w:  and  w2  are  equal  to  zero.  This  implies 
that  the  final  estimate  of  the  position  and  the  initial  estimate 
of  the  position  are  equal.  If  the  three  lines  of  position  were 
obtained  by  sextant  observations,  the  assumed  position  should  be 
the  reference  position.  In  this  case,  the  angles  Oj  ,  ctz  and 
a3  would  be  the  azimuth  angles. 


Appendix  1.  The  Program  COMP 

COMP  is  a  program  with  which  to  implement  the  procedure 
developed  in  Section  II  for  generating  a  composite  estimate  of 


the  position  of  an  object  from  two  or  more  independent  estimates 
of  the  position  of  the  object.  The  independent  estimates  must  be 
specified  in  terms  of  coordinates  in  a  rectangular  coordinate 
system  whose  positive  y-exis  is  north  and  positive  x-axis  is 
east.  In  addition,  a  confidence  region  for  each  estimate  must  be 


specified  in  terms  of  the  direction  S  of  the  major  axis  of  a 
bounding  ellipse  and  either  its  size  k  or  its  containment 


probability  p  (confidence) . 


The  program  listing  follows: 


10  CLS  :  PI  =  4  *  ATN(l) 

20  INPUT  "print  data  yes/no  (y/n)";  A$ 

30  IF  A$  =  "N"  OR  A$  =  "n"  THEN  P$  =  "N":  GOTO  60 
40  IF  A$  =  "Y"  OR  A$  =  "y"  THEN  P$  =  "Y"  ELSE  GOTO  20 

50  LPRINT  :  LPRINT  :  LPRINT  :  LPRINT  "Composite  Position  Estimation  Program:  COMP. BAS” 

60  PRINT  :  INPUT  "estimates  to  be  combined";  NO 
70  IF  P$  =  "N"  THEN  GOTO  90 

80  LPRINT  :  LPRINT  "estimates  to  be  combined  ",  NO 

90  INPUT  "elliptical  confidence  regions  by  size  or  probability  (s/p)";  F$ 

100  IF  F$  =  "S"  OR  FJ  =  "s"  OR  FJ  =  "P"  OR  F$  =  "p"  THEN  GOTO  110  ELSE  GOTO  90 

110  IF  F$  =  "P"  OR  F$  =  "p”  THEN  B$  =  "probability"  ELSE  B$  =  "size" 

120  FOR  10  =  1  TO  NO 

130  PRINT  :  INPUT  "position  estimate  x-coordinate";  X 

140  INPUT  "position  estimate  y-coordinate";  Y 

150  INPUT  "containment  ellipse  major  axis  direction";  DO 

160  IF  DO  >=  180  OR  00  <  0  THEN  GOTO  150 

170  INPUT  "ellipse  major  axis  length";  MA:  SMA  =  MA  /  2 

180  INPUT  "ellipse  minor  axis  length";  MI:  SMI  =  MI  /  2 

190  IF  F$  =  "P"  OR  F$  =  "p"  THEN  GOTO  220 

200  INPUT  "containment  ellipse  size";  SO:  P0  =  1  -  EXP(-S0  *  SO  /  2) 

210  PRINT  "containment  probability  =  ";  P0:  GOTO  240 

220  INPUT  "containment  probability";  P0:  SO  *  SQR ( -2  *  L0G(1  -  P0) } 

230  PRINT  "containment  ellipse  size  =  SO 
240  IF  P$  =  "N"  THEN  GOTO  320 

250  LPRINT  ;  LPRINT  :  LPRINT  "position  estimate  x-coordinate",  X 

260  LPRINT  "position  estimate  y-coordlnate",  Y 

270  LPRINT  "ellipse  major  axis  direction",  DO 

280  LPRINT  "ellipse  major  axis  length  ",  MA 

290  LPRINT  "ellipse  minor  axis  length  ",  MI 

300  LPRINT  "ellipse  size  ”,  SO 

310  LPRINT  "containment  probability  ",  P0 

320  SX  «  SMA  /  SO:  SY  =  SMI  /  SO 

330  DO  =  DO  *  PI  /  180 

340  U  «  SX  *  SX  *  SIN(00)  *  SIN(00)  +  SY  *  SY  *  COS(DO)  *  C0S(D0) 

350  V  *  SX  *  SX  *  C0S(D0)  *  COS(OO)  +  SY  *  SY  *  SIN(DO)  *  SIN(DO) 

360  W  *  (SX  *  SX  -  SY  *  SY)  *  SIN(DO)  *  C0S(D0) 

370  Z=U*V-W*W:  A  =  A  +  V  /  Z:  8  *  B  -  W  /  Z:  C  =  C  +  U  /  Z:  D  =  D  +  V/  Z*X-W/Z*Y 
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380  E  =  E-  W/  Z*X  +  U/  Z*Y:F  =  F  +  V/  Z*V/Z*U:G  =  G  +  W/  Z*W/Z*U 

390  H  =  H-  V/  Z*W/Z*U:  I=I+W/Z*W/Z*V:J=J+U/Z*U/Z*V 

400  K=K-V/Z*U/Z*V:  L  =  L-  V/  Z*W/Z*W 

410  M  =  M+(V/Z*U/Z  +  W/  Z*W/Z)*U:N  =  N-  W/  Z*U/Z*W 

420  NEXT  10 

430  Q=A*C-B*B:X=(C*0-B*E)/Q:Y=(A*E-B*D)/Q 

440  R=(C*C*(F+I+2*L)-2*C*B*(H+K+M)+B*B*(6+J+2*  N))  /  (Q  *  Q) 

450  S  »  (A  *  A  *  (G  +  J  +  2  *  K)  -  2  *  A  *  B  *  CH  +  1C  +  M)  +  6  *  8  *  (F  +  I  +  2  *  L))  /  (Q  *  Q) 

460  T  =  {(A  *  C  +  B  *  B)  *  (H  +  K  +  M)  -  C  *  B  *  (F  +  I  +  2  *  L)  -  B  *  A  *  (G  +  J  +  2  *  N) )  /  (Q  *  Q) 

470  CLS  :  PRINT  "composite  estimate  x-coordinate  =  X 
480  PRINT  "composite  estimate  y-coordinate  =  Y 
490  IF  PS  =  ”N"  THEN  GOTO  520 

500  LPRINT  LPRINT  :  LPRINT  :  LPRINT  "composite  estimate  x-coordinate",  X 
510  LPRINT  "composite  estimate  y-coordinate",  Y 
520  IF  S  =  R  THEN  A  =  0:  GOTO  600:  REM  a  circle 

530  A  =  .5  *  ATN(2  *  T  /  (S  -  R ) ) 

540  R1  =  R  *  COS(A)  *  COS(A)  -  2  *  T  *  COS(A)  *  SIN(A)  +  S  *  SIN{A)  *  SIN(A) 

W  Si  =  R  *  SIN(A)  *  SIN(A)  +  2  *  T  *  COS(A)  *  SIN(A)  +  S  *  COS(A)  *  COS(A) 

560  A  =  A  *  180  /  PI 

570  IF  R1  >  SI  THEN  A  =  A  +  90:  GOTO  590 
580  RT  =  Rl:  R1  =  SI:  SI  =  RT 

590  IF  A  <  0  THEN  A  =  A  +  180 

600  PRINT  :  PRINT  :  INPUT  "confidence  region  by  size  or  probability  (s/p)";  G$ 

610  IF  GS  =  "S"  OR  G$  =  "s"  THEN  GOTO  630 

620  IF  G$  =  "P"  OR  G$  =  "p"  THEN  GOTO  650  ELSE  GOTO  600 

630  PRINT  :  INPUT  "containment  ellipse  size";  SO:  PO  =  1  -  EXP ( -SO  *  SO  /  2) 

640  PRINT  :  PRINT  "containment  probability  =  ";  PO:  GOTO  670 

650  PRINT  :  INPUT  "containment  probability";  PO:  SO  =  SQR(-2  *  LOG(l  -  PO)) 

660  PRINT  :  PRINT  "containment  ellipse  size  =  ";  SO 

670  PRINT  "ellipse  major  axis  length  =  ";  2  *  SO  *  SQR(Rl) 

680  PRINT  "ellipse  major  axis  direction  =  ";  A 

690  PRINT  "ellipse  minor  axis  length  =  ";  2  *  SO  *  SQR(Sl) 

700  IF  P$  =  "N"  THEN  GOTO  760 

710  LPRINT  :  LPRINT  "containment  ellipse  size  ",  SO 

720  LPRINT  "containment  probability  ",  PO 

730  LPRINT  "ellipse  major  axis  length  ",  2  *  SO  *  SQR(Rl) 

740  LPRINT  "ellipse  major  axis  direction  ",  A 

750  LPRINT  "ellipse  minor  axis  length  ",  2  *  SO  *  SQR(Sl) 

760  PRINT  :  PRINT  :  INPUT  "continue  or  quit  ( c/q ) " ;  G$ 

770  IF  GJ  =  ”C"  OR  G|  =  "c"  THEN  GOTO  600 

780  IF  G$  =  "Q"  OR  G$  =  "q"  THEN  GOTO  790  ELSE  GOTO  770 

790  END 
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Appendix  2.  The  Program  PEST 


PEST  is  a  program  with  which  to  implement  the  procedure 
developed  in  Section  III  for  determining  a  position  estimate  for 


an  object  from  a  set  of  observed  bearings  from  or  on  two  or  more 


reference  points  of  known  positions.  In  the  development,  the 
reference  points  are  referred  to  as  stations.  The  bearing  error 


associated  with  an  observation  must  be  specified  for  each 


station.  In  addition,  the  position  of  each  station  must  be 


specified  in  terms  of  its  bearing  and  range  from  a  single 


reference  point  which  may  be  one  of  the  reference  points.  The 


program  generates  an  estimate  for  the  position  of  an  object  and  a 


confidence  region  for  the  estimate  which  is  specified  in  terms  of 


the  direction  S  of  the  major  axis  of  a  bounding  ellipse  and 


either  its  size  k  or  its  containment  probability  p 


(confidence) . 


The  program  listing  follows: 


10  CIS  :  PI  =  4  *  ATN(l) :  DIM  A(7):  I  =  0 
20  INPUT  "print  data  (y/n)”;  AS 
30  IF  A$  =  "N"  OR  A$  =  "n"  THEN  P$  =  "N":  GOTO  60 

40  IF  A$  =  "Y"  OR  AS  =  "y"  THEN  PS  =  "Y"  ELSE  GOTO  20 

50  LPRINT :  LPRINT :  LPRINT  :  LPRINT  "Position  Estimation  Program:  PEST. BAS" 
60  INPUT  "bearings  on  object  or  bearings  on  stations  (o/s)";  BS 
7)  IF  BS  =  "0"  OR  BS  =  "o"  THEN  K  =  0:  GOTO  120 

80  IF  BS  =  "S"  OR  BS  =  "s"  THEN  K  =  1  ELSE  GOTO  60 

90  IF  PS  =  "N"  THEN  GOTO  120 

100  IF  K  =  0  THEN  BS  =  "  object"  ELSE  BS  =  "  stations" 

110  LPRINT  :  LPRINT  "bearings  on"  +  BS:  LPRINT 
120  CLS  :  PRINT  "input  number  ”  +  STRS(I  +  1) 

130  INPUT  "observed  bearing  in  decimal  degrees";  P 

140  INPUT  ’’station  bearing  in  decimal  degrees”;  Q 

150  INPUT  'staiion  range”;  R 

160  INPUT  ’’bearing  error  in  decimal  degrees”;  BE 

170  IF  PS  =  "N”  THEN  GOTO  230 

180  LPRINT  :  LPRINT  "input  number  ’’  +  STRS(I  +  1) 

190  LPRINT  "observed  bearing  in  decimal  degrees".  R 
200  LPRINT  "station  bearing  In  decimal  degrees”,  Q 
210  LPRINT  "station  range  in  chosen  units".  R 
220  LPRINT  "bearing  error  in  decimal  degrees",  BE 
230  IF  K  «  1  THEN  P  *  P  +  180 

240  P  *  P  *  PI  /  180:  Q  =  Q  *  PI  /  180:  BE  *  BE  *  PI  /  180 
250  IF  I  -  2  THEN  GOTO  350 

260  I  ■  I  +  1:  A( I  -  1)  =  P:  A(I  +  1)  =  Q:  A(I  +3)  *  R :  A(I  +  5)  *  BE 
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270  IF  I  =  1  THEN  GOTO  120 

280  X  =  A(4)  *  SIN(A( 2 )  -  A(0) ) :  Y  =  A(5)  *  SIN(A(3)  -  A(1 ) ) :  Z  =  SIN(A(1)  -  A(0)) 

290  IF  Z  =  0  THEN  GOTO  900 

300  U  =  (X  *  S I N ( A ( 1 ) )  -  Y  *  SIN(A(0) ) )  /  Z:  V  =  (X  *  C0S(A(1))  -  Y  *  COS(A(0) ) )  /  Z 
310  FOR  J  =  0  TO  1 

320  P  =  A(J):  Q  =  A ( J  +  2) :  R  =  A ( J  +  4) :  BE  =  A(J  +  6):  GOSUB  940 
330  NEXT  J 
340  GOTO  360 
350  GOSUB  940 

360  CIS  :  INPUT  "continue  data  input  or  generate  an  estimate  (c/g)”;  C$ 

370  IF  C$  =  ”G"  OR  C$  =  "g”  THEN  GOTO  390 

380  IF  CJ  =  "C"  OR  CJ  =  "c"  THEN  GOTO  120  ELSE  GOTO  360 

390  CLS 

400  F  =  (B  *  B  -  A  *  C):  IF  F  =  0  THEN  GOTO  900 

410  XI  =  U  +  (B  *  E  -  C  *  0)  /  F:  Y1  =  V  +  (A  *  E  -  B  *  D)  /  F:  GOSUB  1000 

420  R2  =  Rl:  B2  =  B1  *  180  /  PI:  IF  K  =  1  THEN  B2  =  B2  -  180:  IF  B2  <  0  THEN  B2  =  360  +  B2 

430  T  =  SGN(B)  *  PI  /  4:  IF  A  =  C  THEN  GOTO  450 

440  T  =  .5  *  ATN(2  *  B  /  (A  -  C)) 

450  G  =  (C  *  COS(T)  *  COS(T)  -  2  *  B  *  COS(T)  *  SIN(T)  +  A  *  SIN(T)  *  SIN(T) )  /  (-F):  G  =  SQR(G) 

460  H  =  (C  *  SIN(T)  *  SH.iw  +  2  *  B  *  COS(T)  *  SIN(T)  +  A  *  COS(T)  *  COS(T) )  /  (-F):  H  =  SQR(H) 

470  IF  H  >=  G  GOTO  490 

480  Z  =  H:  H  =  G:  G  =  Z:  T  =  T  +  PI  /  2 

490  CLS 

500  PRINT  "bearing  estimate  =  B2 
510  PRINT  "range  estimate  =  R2:  PRINT 

520  IF  PJ  =  "N"  THEN  GOTO  550 

530  LPRINT  :  LPRINT  :  LPRINT  "bearing  estimate",  B2 
540  LPRINT  "range  estimate",  R2 

550  INPUT  "confidence  region  by  size  or  probability  (s/p)”;  C$ 

560  IF  CJ  *  "S"  OR  CJ  -  "s"  THEN  GOTO  580 

570  IF  C$  =  "P"  OR  C$  =  "p"  THEN  GOTO  650  ELSE  GOTO  550 

580  INPUT  "confidence  region  size";  S 

590  IF  S  <=  0  THEN  GOTO  580 

600  CP  =  1  -  EXP ( -S  *  S  /  2) 

610  PRINT  :  PRINT  "containment  probability  =  CP 

620  IF  P$  =  "N"  THEN  GOTO  720 

630  LPRINT  :  LPRINT  "confidence  region  size",  S 

640  LPRINT  "containment  probability",  CP:  GOTO  720 

650  INPUT  "containment  probability  (confidence)";  CP 

660  IF  CP  >=  1  OR  CP  <=  0  THEN  GOTO  650 

670  S  =  SQR(-2  *  LOG(l  -  CP)) 

680  PRINT  :  PRINT  "confidence  region  size  =  ";  S 

690  IF  PJ  =  "N"  THEN  GOTO  720 

700  LPRINT  :  LPRINT  "containment  probability",  CP 

710  LPRINT  "confidence  region  size",  S 

720  X  =  2  *  S  *  H 

730  PRINT  "major  axis  length  =  ";  X 

740  N  =  T  *  180  /  PI:  IF  N  <  0  THEN  N  =  N  +  180 

750  PRINT  "major  axis  direction  =  ";  N 

760  Y  =  2  *  S  *  G 

770  PRINT  "minor  axis  length  =  ";  Y 

780  PRINT  "containment  ellipse  area  =  ";  PI  *  X  *  Y  /  4 

790  IF  PJ  =  "N"  THEN  GOTO  840 

800  LPRINT  "major  axis  length",  X 

810  LPRINT  "major  axis  direction",  N 

820  LPRINT  "minor  axis  length",  Y 

830  LPRINT  "containment  ellipse  area",  PI  *  X  *  Y  /  4 

840  PRINT  :  PRINT  :  INPUT  "continue  or  quit  (c/q)";  CJ 

850  IF  CJ  =  "C"  OR  C$  «  "c"  THEN  GOTO  870 

860  IF  CJ  =  "Q"  OR  CJ  =  "q"  THEN  END  ELSE  GOTO  840 

870  CLS  :  INPUT  "continue  data  input  or  generate  an  estimate  (c/g)";  CJ 

880  IF  CJ  =  "G"  OR  CJ  =  "g"  THEN  PRINT  :  GOTO  500 

890  IF  CJ  =  "C"  OR  CJ  =  "c”  THEN  GOTO  120  ELSE  GOTO  430 

900  PRINT  "no  solution" 

910  IF  PJ  «  "N"  THEN  GOTO  930 
920  LPRINT  :  LPRINT  "no  solution" 

930  ENO 

940  XI  *  U  -  R  *  SIN(Q) :  Y1  =  V  -  R  *  COS(Q):  GOSUB  1000 
950  W  =  P  -  Bl:  L  -  Rl  *  BE:  IF  L  =  0  THEN  GOTO  900 
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960  6  =  C0S(B1)  /  L:  H  =  SIN(Bl)  /  L 

970  IF  W  >=  PI  THEN  W  =  W  -  2  *  PI:  GOTO  990 

980  IF  U  <=  -PI  THEN  U  =  V  +  2  *  PI 

990  W=W/BE:  A  =  G  *  G  +  A:  8  =  G  *  H  +  B:  C  =  H  *  H  +  C:  D=W*G+D:  E  =  V  *  H  +  E:  RETURN 
1000  R1  =  SQR(X1  *  XI  +  Y1  *  Yl):  IF  R1  =  0  THEN  81  =  0:  RETURN 

1010  IF  ABS (XI  /  Rl)  =  1  THEN  Ml  =  SGN(Xl)  *  PI  /  2  ELSE  Ml  =  ATN(X1  /  R1  /  SQR(1  -  XI  *  XI  /  R1  /  R1 ) ) 

1020  IF  ABS(Y1  /  Rl)  =  1  THEN  81  =  (1  -  SGN(Yl))  *  PI  /  2  ELSE  B1  =  PI  /  2  -  ATN(Y1  /  Rl  /  SQR(1  -  Yl  *  Yl  / 
Rl  /  Rl)) 

1030  IF  Ml  <  0  THEN  81  =  2  *  PI  -  B1 
1040  RETURN 
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